Calcul des profils de
piétons au moyen d’une AFDM et d’une CAH
Chargement des
packages nécessaires
#Chargement des packages
library(FactoMineR) # Pour realiser une analyse multivariee
library(openxlsx) # Pour importer/exporter depuis/vers Excel
library(factoextra) # pour produire des graphiques esthetiques a partir des resultats de l'AFDM
library(ggplot2) # package necessaire au fonctionnement de factoextra
library(sf) # pour manipuler les objets ayant une composante spatiale
library(spatstat) # pour le lissage spatial
library(sp) # pour utiliser la fonction as.owin de spatstat
library(maptools) # pour utiliser la fonction as.owin de spatstat
library(cartography) # pour la representation cartographique et l'habillage des cartes
library(raster) # pour le traitement de donnees matricielles
library(tidyverse) # pour utiliser la fonction "separate"
Import du jeu de
données et résumé statistique
VariablesTypoUltraMarcheurs <- read.xlsx("02_analyse_ultra_marcheurs_v6.xlsx")
rownames(VariablesTypoUltraMarcheurs) <- VariablesTypoUltraMarcheurs$id_hog_pers # pour ajouter l'id des personnes comme entête de ligne
str(VariablesTypoUltraMarcheurs) # Pour connaître le type des variables et avoir un aperçu des valeurs prises pour chaque variable
## 'data.frame': 15482 obs. of 24 variables:
## $ id_hog : chr "P01283" "8856" "5268" "3478" ...
## $ id_hog_pers : chr "P01283-4" "8856-1" "5268-2" "3478-1" ...
## $ ocasional : chr "regular" "regular" "regular" "regular" ...
## $ fq_sem_viaje : num 39765 30081 23469 22180 19883 ...
## $ duracion_prom : num 24.5 36.5 8 23.8 30 ...
## $ calif_exp : num 5 4.81 3.14 5 5 ...
## $ volv_casa : num 19883 15598 10058 11090 9941 ...
## $ trabajo : num 0 0 6705 0 0 ...
## $ escuela : num 0 0 0 0 9941 ...
## $ cuidado : num 19883 10027 3353 11090 0 ...
## $ compras : num 0 4456 3353 0 0 ...
## $ edad : chr "15-39" "15-39" "15-39" "15-39" ...
## $ niv_edu : chr "prim_sec" "univ" "media" "media" ...
## $ genero : chr "muj" "muj" "muj" "muj" ...
## $ difisi : chr "difisi_no" "difisi_no" "difisi_no" "difisi_no" ...
## $ licond : chr "licond_no" "licond_no" "licond_no" "licond_no" ...
## $ estratoCL : chr "s_1-2" "s_1-2" "s_1-2" "s_1-2" ...
## $ tot_pers : num 7 5 4 3 7 5 3 2 5 2 ...
## $ num_moto : num 0 0 1 0 0 1 0 0 0 0 ...
## $ num_carro : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num_bici : num 0 1 1 2 0 0 1 1 0 0 ...
## $ GrandesZones : chr "Sur" "Sur" "Sur" "Sur" ...
## $ f_exp : num 1988 1114 671 584 1988 ...
## $ ZAT_Dest_1er_viaje_desde_casa: num 544 714 351 993 567 663 114 115 547 452 ...
head(VariablesTypoUltraMarcheurs)
## id_hog id_hog_pers ocasional fq_sem_viaje duracion_prom calif_exp
## P01283-4 P01283 P01283-4 regular 39765.36 24.50000 5.000000
## 8856-1 8856 8856-1 regular 30081.17 36.48148 4.814815
## 5268-2 5268 5268-2 regular 23468.56 8.00000 3.142857
## 3478-1 3478 3478-1 regular 22180.43 23.81579 5.000000
## P01283-6 P01283 P01283-6 regular 19882.68 30.00000 5.000000
## 3690-2 3690 3690-2 regular 19044.74 9.00000 5.000000
## volv_casa trabajo escuela cuidado compras edad niv_edu genero
## P01283-4 19882.678 0.000 0.000 19882.678 0.000 15-39 prim_sec muj
## 8856-1 15597.642 0.000 0.000 10027.056 4456.469 15-39 univ muj
## 5268-2 10057.953 6705.302 0.000 3352.651 3352.651 15-39 media muj
## 3478-1 11090.215 0.000 0.000 11090.215 0.000 15-39 media muj
## P01283-6 9941.339 0.000 9941.339 0.000 0.000 6-15 prim_sec muj
## 3690-2 7617.896 0.000 0.000 7617.896 3808.948 40-64 prim_sec muj
## difisi licond estratoCL tot_pers num_moto num_carro num_bici
## P01283-4 difisi_no licond_no s_1-2 7 0 0 0
## 8856-1 difisi_no licond_no s_1-2 5 0 0 1
## 5268-2 difisi_no licond_no s_1-2 4 1 0 1
## 3478-1 difisi_no licond_no s_1-2 3 0 0 2
## P01283-6 difisi_no licond_no s_1-2 7 0 0 0
## 3690-2 difisi_no licond_no s_1-2 5 1 0 0
## GrandesZones f_exp ZAT_Dest_1er_viaje_desde_casa
## P01283-4 Sur 1988.2678 544
## 8856-1 Sur 1114.1173 714
## 5268-2 Sur 670.5302 351
## 3478-1 Sur 583.6956 993
## P01283-6 Sur 1988.2678 567
## 3690-2 Sur 761.7896 663
summary(VariablesTypoUltraMarcheurs) # Résumé des valeurs contenues dans l'objet
## id_hog id_hog_pers ocasional fq_sem_viaje
## Length:15482 Length:15482 Length:15482 Min. : 0.0
## Class :character Class :character Class :character 1st Qu.: 0.0
## Mode :character Mode :character Mode :character Median : 567.6
## Mean : 1069.7
## 3rd Qu.: 1521.6
## Max. :39765.4
##
## duracion_prom calif_exp volv_casa trabajo
## Min. : 0.00 Min. :0.000 Min. : 0.0 Min. : 0.00
## 1st Qu.: 9.00 1st Qu.:3.667 1st Qu.: 93.9 1st Qu.: 0.00
## Median : 15.00 Median :4.857 Median : 295.4 Median : 0.00
## Mean : 21.09 Mean :4.085 Mean : 553.1 Mean : 86.79
## 3rd Qu.: 25.00 3rd Qu.:5.000 3rd Qu.: 756.5 3rd Qu.: 0.00
## Max. :1080.00 Max. :9.000 Max. :19882.7 Max. :7105.89
##
## escuela cuidado compras edad
## Min. : 0.0 Min. : 0.0 Min. : 0.00 Length:15482
## 1st Qu.: 0.0 1st Qu.: 0.0 1st Qu.: 0.00 Class :character
## Median : 0.0 Median : 0.0 Median : 0.00 Mode :character
## Mean : 221.9 Mean : 140.7 Mean : 92.85
## 3rd Qu.: 181.3 3rd Qu.: 0.0 3rd Qu.: 0.00
## Max. :9941.3 Max. :19882.7 Max. :5265.09
##
## niv_edu genero difisi licond
## Length:15482 Length:15482 Length:15482 Length:15482
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## estratoCL tot_pers num_moto num_carro
## Length:15482 Min. : 1.000 Min. :0.0000 Min. :0.0000
## Class :character 1st Qu.: 3.000 1st Qu.:0.0000 1st Qu.:0.0000
## Mode :character Median : 4.000 Median :0.0000 Median :0.0000
## Mean : 3.908 Mean :0.1723 Mean :0.2958
## 3rd Qu.: 5.000 3rd Qu.:0.0000 3rd Qu.:1.0000
## Max. :13.000 Max. :8.0000 Max. :6.0000
##
## num_bici GrandesZones f_exp
## Min. : 0.0000 Length:15482 Min. : 1.00
## 1st Qu.: 0.0000 Class :character 1st Qu.: 49.69
## Median : 1.0000 Mode :character Median : 108.10
## Mean : 0.9474 Mean : 140.62
## 3rd Qu.: 2.0000 3rd Qu.: 189.25
## Max. :50.0000 Max. :1988.27
##
## ZAT_Dest_1er_viaje_desde_casa
## Min. : 0.0
## 1st Qu.: 351.0
## Median : 571.0
## Mean : 615.8
## 3rd Qu.: 794.0
## Max. :1908.0
## NA's :766
Apurement du jeu de
données
VariablesTypoUltraMarcheurs <- na.omit(VariablesTypoUltraMarcheurs) # suppression des individus pour lesquels des valeurs ne sont pas renseignées
VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$num_bici <= 10),] # suppression des individus résidant dans des ménages dans lesquels il y a plus de 10 vélos (situation peu probable, sans doute des erreurs de saisie)
VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$calif_exp <= 5),] # suppression des individus ayant des qualifications d'experiences de trajet superieures à 5
VariablesTypoUltraMarcheurs <- VariablesTypoUltraMarcheurs[which(VariablesTypoUltraMarcheurs$duracion_prom <= 300),] # suppression des individus ayant une durée moyenne de déplacement à pied supérieure à 5 heures (300 minutes)
Suppression des
colonnes inutiles pour l’AFDM
# colonnes d'identification des individus
VariablesTypoUltraMarcheurs$id_hog <- NULL
VariablesTypoUltraMarcheurs$id_hog_pers <- NULL
# colonnes contenant des indications de localisation spatiale
VariablesTypoUltraMarcheurs$ZAT_Dest_1er_viaje_desde_casa <- NULL
# Autres colonnes
VariablesTypoUltraMarcheurs$f_exp <- NULL
View(VariablesTypoUltraMarcheurs) # Resume variables
summary(VariablesTypoUltraMarcheurs)
## ocasional fq_sem_viaje duracion_prom calif_exp
## Length:14218 Min. : 0.0 Min. : 1.00 Min. :1.000
## Class :character 1st Qu.: 0.0 1st Qu.: 10.00 1st Qu.:4.000
## Mode :character Median : 628.5 Median : 15.00 Median :5.000
## Mean : 1109.0 Mean : 20.83 Mean :4.234
## 3rd Qu.: 1548.6 3rd Qu.: 25.00 3rd Qu.:5.000
## Max. :39765.4 Max. :291.25 Max. :5.000
## volv_casa trabajo escuela cuidado
## Min. : 0.0 Min. : 0.00 Min. : 0.0 Min. : 0.0
## 1st Qu.: 111.0 1st Qu.: 0.00 1st Qu.: 0.0 1st Qu.: 0.0
## Median : 323.0 Median : 0.00 Median : 0.0 Median : 0.0
## Mean : 574.5 Mean : 87.98 Mean : 233.2 Mean : 144.2
## 3rd Qu.: 769.5 3rd Qu.: 0.00 3rd Qu.: 223.7 3rd Qu.: 0.0
## Max. :19882.7 Max. :7105.89 Max. :9941.3 Max. :19882.7
## compras edad niv_edu genero
## Min. : 0.00 Length:14218 Length:14218 Length:14218
## 1st Qu.: 0.00 Class :character Class :character Class :character
## Median : 0.00 Mode :character Mode :character Mode :character
## Mean : 96.57
## 3rd Qu.: 3.00
## Max. :5265.09
## difisi licond estratoCL tot_pers
## Length:14218 Length:14218 Length:14218 Min. : 1.000
## Class :character Class :character Class :character 1st Qu.: 3.000
## Mode :character Mode :character Mode :character Median : 4.000
## Mean : 3.886
## 3rd Qu.: 5.000
## Max. :13.000
## num_moto num_carro num_bici GrandesZones
## Min. :0.0000 Min. :0.000 Min. : 0.0000 Length:14218
## 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.: 0.0000 Class :character
## Median :0.0000 Median :0.000 Median : 1.0000 Mode :character
## Mean :0.1678 Mean :0.288 Mean : 0.9254
## 3rd Qu.:0.0000 3rd Qu.:0.000 3rd Qu.: 2.0000
## Max. :8.0000 Max. :5.000 Max. :10.0000
str(VariablesTypoUltraMarcheurs)
## 'data.frame': 14218 obs. of 20 variables:
## $ ocasional : chr "regular" "regular" "regular" "regular" ...
## $ fq_sem_viaje : num 39765 30081 23469 22180 19883 ...
## $ duracion_prom: num 24.5 36.5 8 23.8 30 ...
## $ calif_exp : num 5 4.81 3.14 5 5 ...
## $ volv_casa : num 19883 15598 10058 11090 9941 ...
## $ trabajo : num 0 0 6705 0 0 ...
## $ escuela : num 0 0 0 0 9941 ...
## $ cuidado : num 19883 10027 3353 11090 0 ...
## $ compras : num 0 4456 3353 0 0 ...
## $ edad : chr "15-39" "15-39" "15-39" "15-39" ...
## $ niv_edu : chr "prim_sec" "univ" "media" "media" ...
## $ genero : chr "muj" "muj" "muj" "muj" ...
## $ difisi : chr "difisi_no" "difisi_no" "difisi_no" "difisi_no" ...
## $ licond : chr "licond_no" "licond_no" "licond_no" "licond_no" ...
## $ estratoCL : chr "s_1-2" "s_1-2" "s_1-2" "s_1-2" ...
## $ tot_pers : num 7 5 4 3 7 5 3 2 5 2 ...
## $ num_moto : num 0 0 1 0 0 1 0 0 0 0 ...
## $ num_carro : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num_bici : num 0 1 1 2 0 0 1 1 0 0 ...
## $ GrandesZones : chr "Sur" "Sur" "Sur" "Sur" ...
## - attr(*, "na.action")= 'omit' Named int [1:767] 77 78 185 279 288 539 594 595 596 650 ...
## ..- attr(*, "names")= chr [1:767] "964-1" "964-3" "P00002-2" "4497-1" ...
Calcul de la première
AFDM (Analyse factorielle des Données Mixtes)
res <- FAMD(VariablesTypoUltraMarcheurs, ncp = 10, graph = FALSE)
#Affichage des plans factoriels
plot(res, (choix = c("ind", "var", "quanti", "quali")), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

plot(res, choix = c("var"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

plot(res, choix = c("quanti"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

plot(res, choix = c("quali"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

#Affichage de la fréquence cumulée des valeurs propres
get_eigenvalue(res)
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.408112 10.993911 10.99391
## Dim.2 2.457661 7.927940 18.92185
## Dim.3 1.871238 6.036250 24.95810
## Dim.4 1.661313 5.359074 30.31717
## Dim.5 1.303185 4.203824 34.52100
## Dim.6 1.265837 4.083345 38.60434
## Dim.7 1.120536 3.614631 42.21897
## Dim.8 1.099534 3.546882 45.76586
## Dim.9 1.079988 3.483832 49.24969
## Dim.10 1.040306 3.355826 52.60551
fviz_screeplot(res)

Affichage des
variables pour en placer en supplémentaire
print(as.data.frame(colnames(VariablesTypoUltraMarcheurs)))
## colnames(VariablesTypoUltraMarcheurs)
## 1 ocasional
## 2 fq_sem_viaje
## 3 duracion_prom
## 4 calif_exp
## 5 volv_casa
## 6 trabajo
## 7 escuela
## 8 cuidado
## 9 compras
## 10 edad
## 11 niv_edu
## 12 genero
## 13 difisi
## 14 licond
## 15 estratoCL
## 16 tot_pers
## 17 num_moto
## 18 num_carro
## 19 num_bici
## 20 GrandesZones
Calcul d’une deuxième
AFDM après passage des variables peu discriminantes en
supplémentaire
res <- FAMD(VariablesTypoUltraMarcheurs, ncp = 6, sup.var = c(3, 4, 17), graph = FALSE) ##sup.var signifie qu'on passe les variables calif_exp (3), duracion_prom (4) et num_moto (17) en supplémentaires. ncp = 6 --> on ne garde que 6 axes, qui expliquent déjà 40% de l'inertie totale.
#Affichage des plans factoriels
plot(res, (choix = c("ind", "var", "quanti", "quali")), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes", palette=palette(c("grey","red","blue")), title = "Factor map of individuals (pedestrians) and categories")

plot(res, choix = c("var"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes", title = "Variable's factor map", palette=palette(c("grey","red","blue")))

plot(res, choix = c("quanti"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

plot(res, choix = c("quali"), axes = c(1, 2), lab.var = TRUE, lab.ind = FALSE, autoLab = "yes")

#Affichage de la fréquence cumulée des valeurs propres
get_eigenvalue(res)
## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 3.395012 12.125041 12.12504
## Dim.2 2.452448 8.758744 20.88379
## Dim.3 1.846301 6.593932 27.47772
## Dim.4 1.643386 5.869236 33.34695
## Dim.5 1.281049 4.575174 37.92213
## Dim.6 1.220441 4.358717 42.28084
fviz_screeplot(res)

Calcul de la CAH
(Classification Ascendante Hiérarchique) et graphiques
nbclasses <- 6 # Definition du nombre de classes.
res.hcpc <- HCPC(res, kk=Inf, nb.clust=nbclasses, consol=FALSE, graph = FALSE)
#consol = FALSE signifie qu'on n'applique pas a l'issue de la CAH une consolidation par les k-means. graph = FALSE pour ne pas afficher les sorties graphiques. kk=Inf signifie qu'aucune partition prealable n'est réalisée.
plot(res.hcpc, choice="tree", tree.barplot = FALSE, cex = 0.01)

#Affichage dans le plan factoriel des classes
fviz_cluster(res.hcpc,
shape = 20,
repel = FALSE, # Avoid label overlapping
show.clust.cent = TRUE, # Show cluster centers
palette = "jco", # Color palette see ?ggpubr::ggpar
ggtheme = theme_minimal(),
main = "Factor map",
)

Création d’un tableau
contenant les effectifs par classe et de l’histogramme
correspondant
res.hcpc.dataclust <- as.data.frame(res.hcpc[["data.clust"]]) # conversion du tableau obtenu avec la CAH en dataframe
TableauContingence <- data.frame(res.hcpc.dataclust[,"clust"]) # Préparation d'un tableau des effectifs de chaque classe
names(TableauContingence) <- c("clust")
TableauContingence$Count <- 1
TableauContingence <- aggregate(Count ~ clust, TableauContingence, sum)
View(TableauContingence) # Affichage du tableau d'effectifs de classe
barplot(TableauContingence$Count, names = TableauContingence$clust, main = "Individuals' distribution by cluster") #Affichage des effectifs de chaque classe en histogramme.

Identification des
variables catégorielles ne prenant que deux modalités afin d’alléger la
caractérisation des profils.
# Cette opération est nécessaire pour simplifier ensuite les graphiques
TableauVariablesCategorielles <- Filter(is.character, VariablesTypoUltraMarcheurs)
NbModalitesParVariable <- data.frame()
for (i in colnames(TableauVariablesCategorielles)){
# Recuperation du nombre de modalités pour chaque variable catégorielle
fichier<- length(unique(TableauVariablesCategorielles[[i]]))
NbModalitesParVariable <- rbind(NbModalitesParVariable, fichier)
}
NbModalitesParVariable$Variable <- colnames(TableauVariablesCategorielles)
VariablesDeuxModalites <- NbModalitesParVariable[which(NbModalitesParVariable[,1] == 2),]
VariablesDeuxModalites <- as.character(VariablesDeuxModalites$Variable)
Définition du seuil
de la valeur V.test à retenir pour conserver les variables décrivant les
différentes classes
# Par défaut v.test > |1.96| - Relever le seuil permet de filtrer le nombre de variables qui décriront les classes. On ne garde que celles qui sont les plus significatives
Seuil.v.test <- 3.29
# Avec ce seuil, on accepte un seuil d'erreur à 0,001.
# en savoir plus https://www.medcalc.org/manual/values-of-the-normal-distribution.php
Création des
tableaux décrivant les différentes classes par les variables quanti et
quali et affichage des diagrammes
# On ne garde qu'une modalité pour les variables n'ayant que deux modalités
for (i in 1:nbclasses){
# Récuperation dans des tableaux de la description de chacune des n classes par les variables quantitatives
b <- as.data.frame(res.hcpc$desc.var$quanti[[i]])
b <- b[which(colnames(b) == 'v.test')]
b <- signif(b,3) # pour arrondir les valeurs
# Récuperation dans des tableaux de la description de chacune des n classes par les variables catégorielles
c <- as.data.frame(res.hcpc$desc.var$category[[i]])
d <- as.data.frame(c$label <- rownames(c))
d <- as.data.frame(c$label <- sub("\\=.*", "", c$label))
d <- as.data.frame(c[which(c$label %in% VariablesDeuxModalites), ])
d <- as.data.frame(d[which(d$v.test > 0),])
e <- as.data.frame(c[which(!c$label %in% VariablesDeuxModalites), ])
f <- as.data.frame(rbind(d, e))
g <- as.data.frame(f$label <- NULL)
c <- signif(f,3) # pour arrondir les valeurs
c <- c[which(colnames(c) == 'v.test')]
# Assemblage des tableaux avec les variables quantitatives et qualitatives
h <- as.data.frame(rbind(b,c))
h <- as.data.frame(cbind(h, h$label <- row.names(h)))
colnames(h) <- c("v.test", "label")
h <- as.data.frame(h[order(h$v.test, decreasing = TRUE), ])
h <- as.data.frame(h[which(!h$v.test == "Inf"),])
h <- as.data.frame(h[which(abs(h$v.test) >= Seuil.v.test),])
# Pour tracer l'ensemble des variables significatives sur les diagrammes en barres, il faut au préalable remplacer les V.test infinies (-Inf et Inf). On propose de les remplacer en reprenant les V.test min et max auxquelles on applique un coefficient
j <- h[which(h$v.test != -Inf),] # création tableau sans les valeurs -Inf
min.j <- min(j$v.test) # sélection valeur minimale
h$v.test[which(h$v.test == -Inf)] <- min.j*1.2 # remplacement des -Inf par valeur minimale (multipliée par un coefficient 1.2)
j <- h[which(h$v.test != Inf),] # création tableau sans les valeurs Inf
max.j <- max(j$v.test) # sélection valeur maximale
h$v.test[which(h$v.test == Inf)] <- max.j*1.2 # remplacement des Inf par valeur maximale (multipliée par un coefficient 1.2)
# pour tracer les diagrammes en barres
op <- par(oma=c(0,7,0,0))
barplot((h$v.test), names = row.names(h), border = "white", horiz = TRUE, las = 1, xlim = c(-5-max(abs(h$v.test)), 5 + max(abs(h$v.test))), cex.names=0.8, main = paste("profil",i, sep=""))
}






Récupération de la
variable clust, et réintégration de la variable id_hog_pers au jeu de
données initial
VariablesTypoUltraMarcheurs$profil <- res.hcpc.dataclust$clust
VariablesTypoUltraMarcheurs$id_hog_pers <- row.names(res.hcpc.dataclust)
# utile par la suite pour récupérer les coord. XY du lieu de résidence du ménage
df <- data.frame(x = VariablesTypoUltraMarcheurs$id_hog_pers)
df <- df %>% separate(x, c("id_hog", "id_pers"))
VariablesTypoUltraMarcheurs$id_hog <- df$id_hog
Cartographie des lieux
de résidence des marcheurs suivant leur profil
Importation de la
couche SIG des UTAM (Unités territoriales d’analyse de la mobilité)
UTAM <- st_read("EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3.gpkg")
## Reading layer `EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3' from data source `F:\Rennes\Recherche\Thèses_suivies\Arthur Ducasse\Article_Typo_pietons\EMU2019_UTAM_Bgta_Carto_Variables_Nbre_Completo_v3.gpkg'
## using driver `GPKG'
## Simple feature collection with 132 features and 55 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N
Création d’un semis
de points (objet sf) à partir des coordonnées XY des lieux de résidence
des marcheurs
ResidMarcheurs<- st_as_sf(VariablesTypoUltraMarcheurs, coords=c("X_Hog_Alea","Y_Hog_Alea"), crs = st_crs(UTAM)$srid)
plot(st_geometry(ResidMarcheurs))

Préparation du jeu de
données pour l’élaboration des cartes
# étape nécessaire en vue de créer des cartes lissées donnant à voir la concentration des lieux de résidence des marcheurs suivant leur profil
# pour créer une liste des profils de marcheurs, dans l'ordre croissant des classes
TypoProfilMarcheurs <- VariablesTypoUltraMarcheurs[order(VariablesTypoUltraMarcheurs$profil, decreasing = FALSE), ]
ListProfils <- unique(TypoProfilMarcheurs$profil)
str(TypoProfilMarcheurs)
## 'data.frame': 14218 obs. of 26 variables:
## $ id_hog : chr "10288" "10324" "10496" "10500" ...
## $ ocasional : chr "regular" "regular" "regular" "ocas" ...
## $ fq_sem_viaje : num 2211 257 780 0 0 ...
## $ duracion_prom: num 15.8 10 8 30 45 ...
## $ calif_exp : num 5 4 4 5 5 3 4 5 5 5 ...
## $ volv_casa : num 1005 0 390 130 130 ...
## $ trabajo : num 0 0 0 0 0 ...
## $ escuela : num 0 257 0 0 0 ...
## $ cuidado : num 0 0 0 0 0 0 0 0 0 0 ...
## $ compras : num 1206 0 0 130 130 ...
## $ edad : chr "sup_65" "men_5" "sup_65" "40-64" ...
## $ niv_edu : chr "univ" "prim_sec" "prim_sec" "univ" ...
## $ genero : chr "muj" "hom" "muj" "hom" ...
## $ difisi : chr "difisi_no" "difisi_no" "difisi_si" "difisi_no" ...
## $ licond : chr "licond_no" "licond_no" "licond_no" "licond_si" ...
## $ estratoCL : chr "s_5-6" "s_5-6" "s_5-6" "s_5-6" ...
## $ tot_pers : num 2 4 2 2 2 1 2 2 2 2 ...
## $ num_moto : num 0 0 0 0 0 0 0 0 0 0 ...
## $ num_carro : num 0 1 1 4 4 1 1 1 1 1 ...
## $ num_bici : num 0 1 0 2 2 0 1 0 1 1 ...
## $ GrandesZones : chr "Norte" "Oeste" "Oeste" "Oeste" ...
## $ profil : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ id_hog_pers : chr "10288-1" "10324-4" "10496-1" "10500-1" ...
## $ ID_MZA_EMU : chr "11001191110709" "11001155060105" "11001163030422" "11001163030422" ...
## $ X_Hog_Alea : num 603204 599329 599155 599215 599215 ...
## $ Y_Hog_Alea : num 520164 515043 514540 514599 514599 ...
##Création des objets sf dont seront extrait les toponymes qui seront affichés sur les cartes.
#Pour les lieux de résidence
topores1 <- UTAM[which(UTAM$UTAM_NOMBR == "CHAPINERO" | UTAM$UTAM_NOMBR == "USAQUEN"| UTAM$UTAM_NOMBR == "EL PRADO"),]
topores2 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "CARVAJAL" | UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "CHAPINERO"),]
topores3 <- UTAM[which(UTAM$UTAM_NOMBR == "LA CANDELARIA"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "ALAMOS"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO"| UTAM$UTAM_NOMBR == "CHAPINERO"),]
topores4 <- UTAM[which(UTAM$UTAM_NOMBR == "EL LUCERO"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "20 DE JULIO" | UTAM$UTAM_NOMBR == "CIUDAD USME"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO" | UTAM$UTAM_NOMBR == "CHIA" | UTAM$UTAM_NOMBR == "CAJICA" | UTAM$UTAM_NOMBR == "EL ROSAL" | UTAM$UTAM_NOMBR == "BOJACA" | UTAM$UTAM_NOMBR == "FACATATIVA"| UTAM$UTAM_NOMBR == "LA CALERA" | UTAM$UTAM_NOMBR == "SIBATE"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "SAN CRISTOBAL NORTE"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "SUBA"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "FACATATIVA"| UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "SAN ISIDRO - PATIOS"| UTAM$UTAM_NOMBR == "COTA"| UTAM$UTAM_NOMBR == "TABIO"| UTAM$UTAM_NOMBR == "ZIPAQUIRA"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "LA CANDELARIA"),]
topores5 <- UTAM[which(UTAM$UTAM_NOMBR == "MINUTO DE DIOS"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "CIUDAD USME" | UTAM$UTAM_NOMBR == "20 DE JULIO"| UTAM$UTAM_NOMBR == "TOCANCIPA"| UTAM$UTAM_NOMBR == "GACHANCIPA"| UTAM$UTAM_NOMBR == "LA CALERA"| UTAM$UTAM_NOMBR == "SOPO"| UTAM$UTAM_NOMBR == "TENJO"| UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "MADRID"| UTAM$UTAM_NOMBR == "EL ROSAL"| UTAM$UTAM_NOMBR == "GRAN YOMASA"| UTAM$UTAM_NOMBR == "BOJACA"| UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "TUNJUELITO"| UTAM$UTAM_NOMBR == "AMERICAS"| UTAM$UTAM_NOMBR == "LA CANDELARIA"),]
topores6 <- UTAM[which(UTAM$UTAM_NOMBR == "EL LUCERO"| UTAM$UTAM_NOMBR == "BOSA OCCIDENTAL" | UTAM$UTAM_NOMBR == "SUBA" | UTAM$UTAM_NOMBR == "VENECIA"| UTAM$UTAM_NOMBR == "ENGATIVA" | UTAM$UTAM_NOMBR == "LA URIBE"| UTAM$UTAM_NOMBR == "SOACHA PUEBLO" | UTAM$UTAM_NOMBR == "FUNZA"| UTAM$UTAM_NOMBR == "DIANA TURBAY" | UTAM$UTAM_NOMBR == "GRAN YOMASA"),]
#Sauvegarde des objets sf dans une liste
toporesid <- list(topores1, topores2, topores3, topores4, topores5, topores6)
Production d’une
série de cartes montrant la concentration des lieux de résidence des
marcheurs selon leur profil
#Définition du nombre de classes
nclass <- 5
#création d'un ensemble de palettes de couleurs (on en crée ici plus qu'il n'en faut, pour avoir de la marge si l'on souhaite augmenter le nombre de classes)
pal1 <- hcl.colors(nclass, "YlOrRd", rev = TRUE)
pal2 <- hcl.colors(nclass, "Greens 3", rev = TRUE)
pal3 <- hcl.colors(nclass, "Reds 3", rev = TRUE)
pal4 <- hcl.colors(nclass, "Grays", rev = TRUE)
pal5 <- hcl.colors(nclass, "Oslo", rev = FALSE)
pal6 <- hcl.colors(nclass, "Red-Purple", rev = TRUE)
pal7 <- hcl.colors(nclass, "Inferno", rev = TRUE)
pal8 <- hcl.colors(nclass, "Reds", rev = TRUE)
pla9 <- hcl.colors(nclass, "Blues 3", rev = TRUE)
pal10 <- hcl.colors(nclass, "Purples 3", rev = TRUE)
#Sauvegarde des palettes dans une liste
pal <- list(pal1, pal2, pal3, pal4, pal5, pal6, pal7, pal8, pla9, pal10)
#Définition des limites des UTAM comme emprise pour le lissage
Emprise <- as.owin(as(st_union(st_buffer(UTAM, 1000)), "Spatial"))
#Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(oma=c(3.5,0,3,0)+0.1, mar = c(0, 0.5, 1.2, 0.5))
plot.new()
#Pour découper la fenêtre en 1 ligne et 6 colonnes (6 profils)
par(mfrow=c(3,2))
#Boucle pour produire et cartographier une surface de densité par profil
for (i in ListProfils){
#Récuperation des jeux de données par profil
fichier <- TypoProfilMarcheurs[which(TypoProfilMarcheurs$profil == i),]
#Pour récupérer les coordonnées des lieux de résidence
pts <- st_coordinates(st_geometry(st_as_sf(fichier, coords=c("X_Hog_Alea","Y_Hog_Alea"), crs = st_crs(UTAM)$srid)))
#Pour créer un objet ppp (format spatstat) et y intégrer l'emprise
fichier.ppp <- ppp(pts[,1], pts[,2], window = Emprise)
# pour définir la taille du rayon (en m.)
rayon <- 1000
# pour définir la taille du pixel (en m.)
resol <- 100 # ici 100m x 100m = 1 ha.
# pour calculer la surface de densité (rayon lissage : 1000 m et résolution spatiale de l'image : 100m x 100m = 1 ha)
cartelissee <- density(fichier.ppp, sigma = rayon, kernel = "gaussian", eps = resol)
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
#Découpage du raster sur l'emprise des UTAM
cartelissee.raster <- mask(cartelissee.raster, UTAM)
#Passage des densités à des effectifs (multiplication des densités par la surface du cercle)
values(cartelissee.raster) <- values(cartelissee.raster) * pi*rayon**2
#Définition des seuils de classes (discrétisation en intervalles égaux)
bks <- seq(from = cellStats(cartelissee.raster, stat = min),
to = cellStats(cartelissee.raster, stat = max),
by = (cellStats(cartelissee.raster, stat = max) - cellStats(cartelissee.raster, stat = min)) / nclass)
#Reclassification de la surface lissée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks, include.lowest = FALSE, right = TRUE, dig.lab = 3, ordered_result = FALSE)
#Vectorisation de la surface reclassée (calcul un peu long)
cartelissee.vecteur <- as(rasterToPolygons(cartelissee.reclass, n=4, na.rm=TRUE, digits=12, dissolve=TRUE), "sf")
#Tracer la carte
plot(st_geometry(UTAM), border = "white", bg= "grey90")
typoLayer(
x = cartelissee.vecteur,
var="layer",
col = unlist(pal[as.numeric(i)]),
lwd = 0.1,
border = unlist(pal[as.numeric(i)]),
legend.pos = "n",
add = TRUE)
labelLayer(
x = toporesid[[as.numeric(i)]],
txt = "UTAM_NOMBR",
col= "black",
cex = 0.3,
font = 4,
halo = TRUE,
bg = "white",
r = 0.05,
overlap = FALSE,
show.lines = FALSE)
legendChoro(
pos = "bottomright",
title.txt = "Number of walker residence\nplaces in a 1 km radius.",
title.cex = 0.6,
breaks = bks,
values.rnd = 2,
nodata = FALSE,
col = unlist(pal[as.numeric(i)]),
border = "white",
horiz = FALSE
)
title(main =paste("Cluster",i, sep="-"))
mtext(text = paste0("n = ", nrow(fichier), " walkers"),
side = 3, line = -2, adj = 0.1, cex = 0.7)
}
barscale(
lwd = 1.5,
cex = 0.6,
pos = "bottomleft",
style = "pretty",
unit = "m"
)
north(pos = "topright")
#Pour afficher le titre principal et la source
mtext("Concentration of walkers' home location by profile in the metropolitan area of Bogota", cex = 1.3, side=3,line=1,adj=0.5,outer=TRUE)
mtext(" Sources : Union Temporal STEER & CNC, EMD, 2019 - SDM - radius buffer : 1000 m, resolution : 1 ha", side=1, line=1, adj=0, cex=0.6, outer=TRUE)

Cartographie des lieux
de destination associés au premier déplacement des marcheurs depuis leur
domicile
Importation de la
couche SIG représentant les lieux de première destination (centroïdes
des ZAT - Zones d’analyse du transport)
#Le premier trajet à pied des marcheurs est répété entre 4 et 5 fois dans la semaine, ce qui le rend assez représentatif des déplacements à pied des marcheurs, puisque le deuxième trajet est un peu moins répété, alors que le 3e est bien plus occasionnel.
destino <- st_read("ZAT_PrimerDestino.gpkg")
## Reading layer `ZAT_PrimerDestino' from data source
## `F:\Rennes\Recherche\Thèses_suivies\Arthur Ducasse\Article_Typo_pietons\ZAT_PrimerDestino.gpkg'
## using driver `GPKG'
## Simple feature collection with 14655 features and 31 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 545952 ymin: 453808 xmax: 652779.8 ymax: 600655.7
## Projected CRS: WGS 84 / UTM zone 18N
plot(st_geometry(destino), pch = "+", cex = 0.5)

Création d’une liste
de marcheurs dont la ZAT de destination du premier déplacement est
renseignée
listdestino <- as.character(destino$id_concat)
Sélection des
marcheurs ayant une ZAT de destination connue pour leur premier
déplacement
TypoProfilMarcheurs <- VariablesTypoUltraMarcheurs[VariablesTypoUltraMarcheurs$id_hog_pers %in% listdestino, ]
Jointure pour
récupérer le profil (clust) sur le lieu de destination
destino <- merge(destino, TypoProfilMarcheurs[c(22:23)], by.x = "id_concat", by.y = "id_hog_pers")
Production d’une
deuxième série de cartes représentant la concentration des lieux
associés à la première destination des marcheurs
#Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(oma=c(3.5,0,3,0)+0.1, mar = c(0, 0.5, 1.2, 0.5))
plot.new()
#Pour découper la fenêtre en 2 lignes et 2 colonnes (4 profils)
par(mfrow=c(3,2))
#Définition du nombre de classes
nclass <- 5
#Définition d'une boucle pour produire et cartographier une surface de densité par profil
for (i in ListProfils){
#Récuperation des jeux de données par profil
fichier <- destino[which(destino$profil == i),]
#Récupération des coordonnées des lieux de résidence
pts <- st_coordinates(st_geometry(fichier))
#Création d'un objet ppp (format spatstat) et intégration de l'emprise
fichier.ppp <- ppp(pts[,1], pts[,2], window = Emprise)
#Définition de la taille du rayon (en m.)
rayon <- 1000
#Définition de la taille du pixel (en m.)
resol <- 100 # ici 100m x 100m = 1 ha.
#Calcul de la surface de densité (rayon lissage : 1000 m et résolution spatiale de l'image : 100m x 100m = 1 ha)
cartelissee <- density(fichier.ppp, sigma = rayon, kernel = "gaussian", eps = resol)
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid
#Découpage du raster sur l'emprise des UTAM
cartelissee.raster <- mask(cartelissee.raster, UTAM)
#Passage des densités à des effectifs (multiplication des densités par la surface du cercle)
values(cartelissee.raster) <- values(cartelissee.raster) * pi*rayon**2
#Définition des seuils de classes (discrétisation en intervalles égaux)
bks <- seq(from = cellStats(cartelissee.raster, stat = min),
to = cellStats(cartelissee.raster, stat = max),
by = (cellStats(cartelissee.raster, stat = max) - cellStats(cartelissee.raster, stat = min)) / nclass)
#Reclassification de la surface lissée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks, include.lowest = FALSE, right = TRUE, dig.lab = 3, ordered_result = FALSE)
#Vectorisation de la surface reclassée (calcul un peu long)
cartelissee.vecteur <- as(rasterToPolygons(cartelissee.reclass, n=4, na.rm=TRUE, digits=12, dissolve=TRUE), "sf")
#Tracé de la carte
plot(st_geometry(UTAM), border = "white", bg= "grey90")
typoLayer(
x = cartelissee.vecteur,
var="layer",
col = unlist(pal[as.numeric(i)]),
lwd = 0.1,
border = unlist(pal[as.numeric(i)]),
legend.pos = "n",
add = TRUE)
labelLayer(
x = topodest[[as.numeric(i)]],
txt = "UTAM_NOMBR",
col= "black",
cex = 0.3,
font = 4,
halo = TRUE,
bg = "white",
r = 0.05,
overlap = FALSE,
show.lines = FALSE)
legendChoro(
pos = "bottomright",
title.txt = "Number of destination places\nin a 1 km radius.",
title.cex = 0.6,
breaks = bks,
values.rnd = 2,
nodata = FALSE,
col = unlist(pal[as.numeric(i)]),
border = "white",
horiz = FALSE
)
title(main =paste("Cluster",i, sep="-"))
mtext(text = paste0("n = ", nrow(fichier), " walkers"),
side = 3, line = -2, adj = 0.1, cex = 0.7)
}
barscale(
lwd = 1.5,
cex = 0.6,
pos = "bottomleft",
style = "pretty",
unit = "m"
)
north(pos = "topright")
# Pour afficher le titre principal et la source
mtext("Concentration of walkers' first destination places in the metropolitan area of Bogota", cex = 1.3, side=3,line=1,adj=0.5,outer=TRUE)
mtext(" Sources : Union Temporal STEER & CNC, EMD, 2019 - SDM - radius buffer : 1000 m, resolution : 1 ha", side=1, line=1, adj=0, cex=0.6, outer=TRUE)
